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Abstract. In this work, we investigate the interaction of free and porous 
media flow by large scale lattice Boltzmann simulations. We study the 
transport phenomena at the porous interface on multiple scales, i.e., 
we consider both, computationally generated pore-scale geometries and 
homogenized models at a macroscopic scale. The pore-scale results are 
compared to those obtained by using different transmission models. Two- 
domain approaches with sharp interface conditions, e.g., of Beavers- 
Joseph-Saffman type, as well as a single-domain approach with a porosity 
depending viscosity are taken into account. For the pore-scale simula¬ 
tions, we use a highly scalable communication-reducing scheme with a 
robust second order boundary handling. We comment on computational 
aspects of the pore-scale simulation and on how to generate pore-scale ge¬ 
ometries. The two-domain approaches depend sensitively on the choice of 
the exact position of the interface, whereas a well-designed single-domain 
approach can significantly better recover the averaged pore-scale results. 

Keywords: lattice Boltzmann method; pore-scale simulation; two-domain 
approach; Darcy-Navier-Stokes coupling; interface conditions 


1 Introduction 

Transport phenomena in porous materials are important in many scientific and 
engineering applications such as catalysis, hydrology, tissue engineering and en¬ 
hanced oil recovery. In the past several decades, flow in porous media has been 
studied extensively both experimentally and theoretically. We refer the inter¬ 
ested reader to the textbook J] and the references therein. In porous media 
flow, we usually distinguish between three scales: the pore-scale, the representa¬ 
tive elementary volume (REV) scale and the domain scale. The REV is defined 
as the minimum element for which macroscopic characteristics of a porous flow 
can be observed. Because experimental setups for many practical questions may 
be too expensive or even impossible to realize, numerical simulation of porous 
media flow can be a useful complementary method to conventional experiments. 
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To describe the flow in the bulk of the porous medium, Darcy’s law is com¬ 
monly used in the form 

/iK -1 u = F - Vp, (1) 

where p is the dynamic viscosity of the fluid, K is the permeability tensor of 
the porous medium, F is the body force, and u and p are averaged velocity and 
pressure quantities, respectively. However, when a porous medium and a free 
flow domain co-exist, e.g., in a river bed, there is no uniquely accepted model for 
the transition between the Darcy model and the free flow. Different approaches 
based on two-domain or on single-domain models are available. Using a single¬ 
domain in combination with the Brinkman equation that modifies Darcy’s law 
by a viscous term 

— /r e ffV 2 u + /rK _1 u = F — Vp, (Br) 

allows to model a smooth transition (see e.g. mmy Here p e g is an effective 
dynamic viscosity in the porous region. However, determining appropriate vis¬ 
cosity parameters for the Brinkman model in the transient region is challenging 
|4|5|6j. Furthermore, the penetration of the flow into the porous medium is found 
to depend on the roughness coefficient of the surface; see e.g. iZiaialial- 

Alternatively, one can use a two-domain approach in combination with a 
sharp transmission condition. Considering the (Navier-)Stokes equation in the 
free flow region and the Brinkman (or Darcy) equation in the porous region, the 
interface plays an important role. Proceeding from the experimental investigation 
of Poiseuille flow over a porous medium, Beavers and Joseph m introduced an 
empirical approach that agreed well with their experiment; see also [3]: They 
suggested to use a slip-flow condition at the interface, i.e., the velocity gradient 
on the fluid side of the interface is proportional to the slip velocity. For simplicity, 
we consider a domain for which the interface is aligned with the flow direction. 
The Beavers-Joseph relation is formulated as 


dU 


dz 


z=0+ 



U m ), 


(BJ) 


where z denotes the coordinate perpendicular to the interface, U = U(z) is the 
mean velocity in flow direction, U s is the slip velocity at the interface 2 = 0 + , 
U m is the seepage velocity that is evaluated far from the plane z = 0 in the 
porous region, k is the permeability, and a is a phenomenological dimensionless 
parameter, only depending on the porous media properties that characterize the 
structure of the permeable material within the boundary region which typically 
varies between 0.01 and 5 mm ■ We refer to mm and the references therein 
for the interface coupling of two-phase compositional porous-media flow and 
one-phase compositional free flow. 

In 1971, Saffman [TB] found that the tangential interface velocity is propor¬ 
tional to the shear stress. He proposed a modification of the BJ condition as 


= (U a -U m ) + 0{k). 


Vk_ dU_ 
a dz 


z= 0+ 


(BJS) 
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More than two decades later, Ochoa-Tapia and Whitaker m proposed an al¬ 
ternative modification of the BJ condition which includes the velocity gradient 
on both sides of the interface as 


/^eff 


dJJ_ 

dz 


2 = 0 “ 


dU 

dz 


2 = 0 + 




(OTW) 


Here the jump-coefficient /3 is a free fitting parameter that needs to be deter¬ 
mined experimentally 1181 . Different expressions for the effective viscosity fi e g 
can be found in the literature. For instance, Lundgren [IS] suggested a relation 
of the form fi e g = fi/e. 

All of the interface conditions mentioned above require the a priory knowl¬ 
edge of the exact position of the interface 120121122 ] , which is for realistic porous 
geometries often not the case. Additionally both, single-domain and two-domain, 
homogenized models rely on assumptions whose validity is not automatically 
guaranteed and depend on additional parameters. Traditional experiments to 
validate and calibrate such models are often costly, time consuming and diffi¬ 
cult to set up. On the other hand, modern high performance computers enable 
the development of increasingly complex and accurate computational models 
resolving pore-scale features. Designing highly efficient solvers for partial differ¬ 
ential equations is one of the challenges of extreme scale computing. While finite 
volume/element/difference schemes give rise to huge algebraic systems, lattice 
Boltzmann methods are intrinsically parallel and extremely attractive from the 
computational complexity point of view. Thus fully resolved direct numerical 
simulation based on first principles modeling is not only feasible nowadays but 
also provides an attractive possibility for validation and calibration. As a step in 
this direction, we here carry out a pore-scale simulation of free flow over a porous 
medium. The model porous media geometry is constructed by generating a ran¬ 
dom sphere-packing using an in-house multi-body simulation framework called 
PE [23j. In the pore geometries constructed such, the flow equations are solved 
with full geometrical resolution. This naturally leads to high computational cost 
requiring the use of high end parallel computing. As we will show by perfor¬ 
mance analysis, the in-house lattice Boltzmann solver WALBerla [Mj exhibits 
excellent performance and parallel scalability for these pore-scale simulations. 

We use the results of the direct numerical simulation of flow over and through 
the porous media as reference solution and evaluate several sharp-interface condi¬ 
tions. As a further example, we also use a homogenized lattice Boltzmann model 
as a REV scale simulation and show the capability of this model to reproduce 
the pore-scale results with high accuracy. 


2 Numerical method 

The lattice Boltzmann method (LBM) has been successfully applied to simulate 
porous media flow. The kinetic nature of the LBM enables it for fluid systems in¬ 
volving microscopic interactions, e.g., flow through porous media. Furthermore, 
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its computational simplicity, its amenability to simple and efficient implemen¬ 
tation and parallelization, and its ability of handling geometrically complex do¬ 
mains make it an applicable tool to simulate porous media How on the pore-scale 
1251261271281 . 

The LBM can also be applied to model the fluid flow in porous media at 
the REV scale. The most commonly used models are the Darcy, the Brinkman- 
extended Darcy and the Forchheimer-extended Darcy models. This last approach 
accounts for the flow resistance in the standard LBM by modifying body-force 
or equilibrium terms, leading to the recovery of either Darcy-Brinkmans equa¬ 
tions or generalized Navier-Stokes equations 12TO31'] . The general model of 
porous media flow should consider the fluid forces and the solid drag force in the 
momentum equation |32| . Guo and Zhao [33j proposed a model to include the 
porosity into the equilibrium distribution and added a force term to the evolu¬ 
tion equation to account for drag forces of the medium. The non-linear inertial 
term is not included in the Brinkman model either, and thus, is only suitable 
for low-speed flow. In this approach, the detailed structure of the medium is 
ignored, and the statistical properties of the medium are included directly. 

2.1 The lattice Boltzmann equation 

The LBM originates from the lattice-gas automata method and can also be 
viewed as a special discrete scheme for the Boltzmann equation with discrete 
velocities 

f (x + e k At, t + At) - f (x, t) = fi(x, t) + F k At, (2) 

where e k is the particle velocity and fl(x, t) is the collision operator. For the 
three dimensional lattice model D 3 Q W , f(x, t) = (/o(x, t), /i(x, t), ..., /is(x, t)) T 
is a 19 dimensional vector of distribution function. Fk is the force that acts as 
a source term to drive the flow. 

It is very common to use the Bhatnagar-Gross-Krook (BGK)[33] model that 
features a single-relaxation-time (SRT) approximation for the collision operator. 
However, it has been shown that using the SRT leads to nonphysical viscosity 
dependence of boundary locations and also suffers from poor stability properties 
pm| . Here, we use the TRT collision operator in which the relaxation time of 
the symmetric and anti-symmetric components of the distribution function are 
separated. For an in-depth discussion of the TRT model we refer to [37138139] . 
As proposed by Ginzburg EZj, the TRT model uses two relaxation rates ui + and 
u>~ where is used for even order moments, and uj~ is used for odd order 
moments 

f2(x,f) = -u + (f + (x,f) - f e<? ’+(x,f)) -u~ (f"(x,t) - f e « “(x,t)) , (3) 

and 

r- j- f k f k r — f k f k / . \ 

Jk 2 ’ ^ k ~ 2 ' V*' 

Here k is the index of the discrete velocity opposite to the one associated with the 
index k. The first eigenvalue is set to l/w + = 3^+ 0.5 and the second eigenvalue 
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ui is a free parameter. Due to stability reasons, oj has to be selected in (0, 2). 
The equilibrium distribution function f eq (x, t n ) for incompressible flow is given 
by m 

/® q (x, t n ) = w k {Sp + p 0 [c~ 2 e k ■ u + \c~ A {e k ■ u) 2 - \c~ 2 u • u] } , (5) 

where w k is a set of weights normalized to unity, p = po + Sp. Here Sp is the 
density fluctuation, and po is the mean density which we set to po = 1 . c s = 
Ax/ (x/3 At) is the lattice speed of sound, while Ax denotes the lattice cell width. 
The macroscopic values of density p and velocity u can be calculated from f as 
zeroth and first order moments with respect to the particle velocity, i.e., 

^=E t=0 ^ u = ^~ E fc=0 e ^- 

In a lattice Boltzmann scheme, we typically split the computation into a 
collision and a streaming step that are given as 


fk (x, t n ) - /fc(x, t n ) = fi(x, t) + F k At, (collision) 

/fc(x + e k At, t n+ 1 ) = /fc(x, tn), (streaming) 

respectively, for k = 0,..., 18. The execution order of these two steps is arbitrary 
and may vary from code to code for implementation reasons. 

In addition, for linear steady flow, it has been demonstrated [35] that most of 
the macroscopic errors/quantities of the TRT depend on A = (Af — \) — 5 ) 

the so-called magic parameter that includes the spatial error, stability, best 
advection and diffusion. The choice A = j is suggested as a suitable value 
for porous media simulations. Another choice, namely A = -A, yields the exact 
location of bounce-back walls in case of Poiseuille flow in a straight channel 

mm . 

2.2 Boundary conditions 

In this study, two types of boundary conditions are used for the pore-scale sim¬ 
ulation. The first one is a no-slip wall condition and the second one is a periodic 
pressure forcing that is applied to drive the flow by a pressure gradient. The sim¬ 
plest scheme to imply no-slip boundary conditions in lattice Boltzmann is the 
simple bounce-back (SBB) operator. In this scheme, the wall location is repre¬ 
sented by a staircase approximation, and the no-slip boundary is satisfied by the 
bounce-back phenomenon of a particle reflecting its momentum upon collision 
with a wall. Hence, the unknown distribution function is calculated as: 

fki x fi>tn+l) = fk(xfi,t n ). (7) 

where we take the values f k after collision but before streaming on the right 
hand side. However, the staircase approximation is not appropriate for complex 
geometries where more accurate results are required even for a low resolution of 




6 


Lecture Notes in Computer Science: Authors’ Instructions 


the boundary. Hence, the central linear interpolation (CLI) scheme which yields 
a higher accuracy at moderately increased computational cost is our preferred 
choice. 

In the CLI scheme [38] three particle distribution functions are needed at 
two fluid nodes adjacent to the solid node, i.e., 

fk( X fl I tn+l) = T+2§ fk( X f2 > tn) ~ T+§f fk( x fl , tn) "b fk{xf 1: t n ). (8) 

while q = | Xf 1 — x w \/\Xf 1 — Xb\ defines a normalized distance of the first fluid 
node to the wall, xf 1 and Xf 2 are the first and second fluid neighbor cells in the 
direction of k, respectively. We use the A = ^ which the CLI has been shown 
to have second order accuracy m- 

3 Large scale simulations 

In this study, we use the WALBerla software framework [24l42j that provides 
a highly optimized implementation of the TRT model that is about as fast as 
the SRT model. We refer to 03 ], where scalability of WALBerla to more than 
10 12 lattice cells and almost 500000 cores has been demonstrated. Here, different 
from the SBB boundary condition, we use the CLI scheme that must access two 
neighboring fluid cells. In WALBerla, this situation is handled by extra ghost- 
layer exchanges, i.e., by communicating an extended set of distribution functions 
to neighboring processors. This results in an additional communication in case 
of massively parallel simulation runs. 

To demonstrate the parallel scalability and efficiency of the WALBerla 
framework in the context of a porous media simulation, we first perform a weak- 
scaling study. Here we use a lattice of 151 3 cells per core and embed into this 
grid a sphere with a diameter of 90 cells. The results have been obtained on the 
LIMA cluster at RRZE0 which has 500 compute nodes. Each node consists of 
two Intel Xeon 5650 ” Westmere” chips so that each node has 12 cores running at 
2.66 GHz. We conduct scalability tests ranging from one node to 64 nodes. This 
setup results in 2.64 x 10 9 cells for the largest run including 768 spherical obsta¬ 
cles. Fig. [D displays the weak-scaling results using the TRT kernel. Fig. |T] shows 
the mega lattice updates per second (MLUPS) for the SBB and CLI boundary 
schemes. The results do not only confirm that the code scales very well, but 
also that the MLUPS value per core compares favorably with many other LBM 
implementations [44145146!] . 

We point out that achieving a good scaling behavior becomes more challeng¬ 
ing when the node performance is already high, but that a high performance 
on each node is a fundamental (though sometimes neglected) prerequisite for 
achieving good overall performance. Thanks to both, the meticulously optimized 
WALBerla kernels on each node, combined with the carefully designed commu¬ 
nication routines, the MLUPS value per core is high and stays nearly constant 
while the number of cores is increased. Note that the CLI boundary condition 

1 https: //www.rrze.fau.de / dienste/arbeiten-rechnen/hpc/systeme 









LBMs for free and porous media flow 


7 


causes a slowdown of about 10% in comparison to the SBB boundary condition, 
which is the fastest scheme. The slowdown of the performance while using the 
CLI is due to the additional time that is needed for the communication and 
the higher complexity of the boundary condition compared to SBB. However 
the higher accuracy of the CLI m compared to the SBB allows in complex 
application to use a coarser resolution of the simulation domain. 



Fig. 1. Measured MLUPS per core for Weak scaling on LIMA-Cluster using 151 3 cells 
per core. 


3.1 Pore-scale simulation with a porous medium generated by a 
particle simulation 

To construct a porous structure, we use the in-house multi-body dynamics frame¬ 
work PE m, The PE can simulate the motion of rigid bodies and their interac¬ 
tion by frictional collisions. Here we use this functionality to generate a random 
sphere packing by letting random spheres fall into the simulation domain from 
the top. After the spheres have come to a rest, their position is fixed and their 
geometry defines the solid matrix of a porous structure. The pore space is then 
resolved by a lattice Boltzmann grid. 

The particles have possibly different radii that vary up to 50% of a mean 
diameter and for each sphere it is chosen randomly. For the fluid flow simulation 
using the LBM, the TRT collision operator and the CLI solid boundary condition 
are used. This combination is fast, has second order accuracy and shows no 
viscosity-dependency. 

First, we test the influence of the cell size on the averaged stream-wise ve¬ 
locity. To do so we increase the diameter D of the spheres from 4 to 48 and keep 
R eD = Um “* D constant. The domain has two walls at the top and bottom, and 
periodic boundary conditions are applied at stream-wise and span-wise direc¬ 
tions. A constant pressure drop drives the flow, and the data are set such that 
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Ren — 2. The simulation result is presented as a planar average of stream-wise 
velocity in Fig. [2] while it is normalized based on the maximum velocity and the 
height of the channel. As it can be seen, the results converge and for D = 32 lat¬ 
tice units, and a further increase of the resolution does not significantly change 
the results. It is worth to note that in the porous region a coarse lattice can be 
used and that only the transient region requires a better resolution. 



Fig. 2. Planar average stream-wise velocity for different grid sizes, ReD — 2. 


Fig. [3(b)] shows the planar average stream-wise velocity for different Re num¬ 
bers. To change the Re number, the viscosity and particles diameter are kept 
constant while the pressure gradient is changed to adjust the flow velocity. The 
results show that for slow flow, the velocity in the porous region is considerably 
higher than for fast flow. When the Re number of the flow increases, the position 
of the maximum velocity shifts toward to the top wall. This observation results 
from a boundary layer effect; when the flow velocity is high in the free flow, the 
penetration to the porous region is less, therefore, the position of the maximum 
velocity changes. 


In Fig. 3(b)| we observe a small deviation in the velocity profile close to the 
bottom wall in the porous region. This is because of the high porosity close to 
the wall, where spherical particles are on a flat plane, see Fig. 3(a) Consequently 
a higher permeability can be found in this region, and the flow will accelerate 
because the resistance against the pressure difference is lower than in the interior 
of the porous medium. Therefore, to evaluate the existing models without this 
effect and having a more uniform porosity in the porous region, a different set-up 
structure is chosen. The bottom plate of the particle simulation is placed about 
one particle size below the bottom wall of the fluid flow simulation. With this 
structure the porosity does not have the effect of placing a sphere on the wall, 
and therefore we have approximately a uniform permeability distribution in the 
porous media. 
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(a) pore geometry 


Porosity 

0 0.2 0.4 0.6 0.8 1 



(b) planar average of stream-wise velocity 


Fig. 3. Flow over mono-sized particles for different Re numbers. 


The results of this pore-scale simulation are taken as a reference solution. 
Here, we use 1274 particles with diameters in the range of 16-48 cells. The flow 
is driven by a pressure difference of 10“ 6 (in lattice units), and the simulation is 
run until the flow reaches the steady state. The planar average of the stream-wise 
velocity is depicted in Fig. [I] 



Porosity 



(b) planar average of stream-wise velocity 


Fig. 4. pore-scale simulation of free flow over porous media. 
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3.2 Evaluation of different interface conditions 

In this subsection, we evaluate different two-domain approaches. All interface 
conditions under consideration have parameters for which no explicit relation 
is known. In the BJ and BJS models, the slip coefficient, a, is unknown, while 
in the OTW model, the jump coefficient /3 and the effective viscosity /z e ff are 
unknown and in the Br model, the effective viscosity fi e ft is unknown. 

By using the DNS solution, we calculate the optimal value for the unknown 
parameters. The domain that is used is a channel which is periodic in stream-wise 
and span-wise directions (Fig. [5]). A free fluid flows on the top of a porous media. 
To have a good comparison, all of the flow properties are non-dimensionalized. 


z* 



6 pen 


6 porous 


Fig. 5. Schematic of the simulation domain and averaged velocity profile in the open 
and porous regions. 


The value of the interface velocity Umti can be directly obtained from the 
averaged velocity profile of the DNS. In order to obtain the velocity gradient 
on the open and porous sides, curve fitting techniques are used to approximate 
the velocity profile close to the interface. The velocity profile on the open side 
can be well approximated by a polynomial curve and on the porous side, the 
velocity profile can be approximated by an exponential curve. Permeability and 
seepage velocity (Darcy velocity) can be calculated from the velocity profile far 
from the interface in the porous medium. Given this, the unknown variables can 
be calculated from the Eqs. amp , am and (lOTWIl . However, to do so, the 
exact position of the interface should be defined which in real applications is 
nearly impossible. 

To find out how the additional parameters of the interface conditions affect 
the results, a two-donrain approach is chosen and solved analytically. For the 
free flow region, the Stokes equation is used and for the porous region, the 
Brinkman’s equation is chosen. The permeability is calculated from the DNS 
result far enough from the interface inside the porous region. In Fig. [Gj we 
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depict the planar average stream-wise velocity which is normalized based on the 
maximum velocity in the DNS solution. 






Fig. 6. Analytical solution for the velocity profile, which is normalized by the maximum 
velocity of the DNS solution, by different interface models. 


As it can be seen in Fig. 6(a) 

Meff 


in the Brinkman model by increasing the 
viscosity ratio, J = the maximum velocity decreases and produces a dis¬ 
continuity in the shear stress over the interface. In the OTW model (Fig. |6(b~)j ), 
negative values of /3 do not influence the result significantly, however, positive 
value of /3 have a strong impact on the maximum velocity as well as on the slip 
velocity on the interface. Fig. |6(c)j and Fig. |6(d)] slrow the results for the BJ and 
the BJS interface conditions. It can be observed that there is almost no differ¬ 
ence between these two models for low Re number flows. In both these cases, 
the maximum velocity decreases if a increases. A small value of a results in a 
considerably larger maximal velocity than in the two other cases. 

Quite often two-domain models result in discontinuities in the stress at the 
interface. Thus the a priori knowledge of the position of the interface is crucial. 
One possibility to fix the position of the interface is to take the location where 
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the porosity reaches the limit value one, i.e., y = 0.756. However fitting of 
the DNS velocity profile shows that only up to y = 0.722, the curve is fitted 
well by an exponential function. More precisely, y = 0.48423 • exp(0.31195 x) — 
0.48236 • exp(0.3131 x) yields a root mean squared error of 5.736 ■ 10~ 6 . The pure 
fluid flow velocity profile is fitted to a 2nd order polynomial resulting in y = 
(1.9593e — 3) + (2.78421e — A)x— (4.48066e — 6)x 2 with a root mean squared error 
of 9.5815-10” 6 . This observation motivates an alternative choice of the interface 
position. Calculating the slip coefficient and the jump coefficient for these two 
positions, we find for y = 0.756, a = 0.3163, /3 = —2.8397 and for y = 0.722, 
a = 0.31645 and /3 = —2.8397. However, as it can be seen in Fig. [3 even with 
the parameters which are extracted from the DNS results, the considered two- 
domain approaches cannot represent accurately the DNS solution. Comparing 
Figs. |7(a)| and |7(b)| shows that the two-domain approaches depend strongly on 
the interface position and more sophisticated criteria for defining the interface 
location are required to obtain better matching results. 




(a) Interface at y=0.756 (b) Interface at y=0.722 

Fig. 7. Normalized velocity profile of the one-domain approaches in compare to the 
DNS solution; a) interface at y=0.756, b) interface at y=0.722. 


4 Comparsion of a homogenized LBM with the pore-scale 
LB simulation 

Different models for isothermal incompressible fluid flow in porous media are 
proposed by several groups. In this work, we use the generalized lattice Boltz¬ 
mann model (GLBM) for porous media introduced in [33], which is applicable 
for a medium with both a constant and a variable porosity. The model can be 
expressed by the following generalized Navier-Stokes equation: 

V ■ u = 0 

+ (u ■ V) 0) = --S7 (ep) + i/ eff V 2 u + F, 


du 

~dt 


(9) 

( 10 ) 
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where p is the fluid density, u and p are the volume-averaged velocity and pres¬ 
sure, respectively, u e s is the effective viscosity, and e is the porosity. The total 
body force F caused by the presence of a porous medium and other external 
force fields is given by 


F = -^u- 7 =|u|u + £ G, (11) 

where v is the shear viscosity of the fluid that is not necessarily the same as 
r'eff, G is the body force induced by an external force, cp is the Forchheimer 
coefficient that depends on the porous structure, and K is permeability of the 
porous media. The first and the second terms on the right hand side of Eq. 
m are the linear Darcy and non-linear Forchheimer drags due to the porous 
medium, respectively. The quadratic nature of the non-linear resistance makes 
it negligible for low-speed flows, but is more noteworthy in hindering the fluid 
motion for high-speed flows, i.e., high Re number and high Da number flows. 

Firstly to validate the generalized model for flow over a porous medium, 
we choose a simple Couette flow. The lower-half of the channel of width H is 
filled with a porous medium with a porosity of e, the stream-wise and span- 
wise boundaries are periodic, and the top wall of the channel is moving with a 
constant velocity of uq. Then, the steady state velocity in this channel satisfies 
the equation 


ev 


CCf 


Z'effV 2 !! - ^=|u|u + eG = 0, 

R V K 

while the walls of the channel are modeled by a no-slip condition. 


( 12 ) 




Fig. 8 . Velocity profile of the Couette flow for different viscosity ratios J = p e /p, in 
comparison with the approximate analytical solution of Eq. G3, a) global system; b) 
zoom into the region near the interface 


Fig. [5] shows the velocity profile for the Couette flow with different viscosity 
ratios J (= p e /p) and compared to a semi-analytical solution for Re = 0.1 and 
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Da = 0.00012. In the Stokes regime for low Da number, [TB] reported that the 
velocity profile in the free flow is linear and exponentially decaying in the porous 
region. More precisely the semi-analytic solution can be written as: 


u x (y) 


rKa + ea(y- H/2) H/2 <y<H 

rKae^y- 11 / 2 ') 0 < y < H/2 


(13) 


while 


2 u 0 


2rK + eH ’ y/i'esk ’ 

and uq is the lid’s velocity. The simulation result shows excellent agreement with 
the analytical solution for both viscosity ratios. 

Secondly, we apply the generalized model to a problem with no sharp interface 
and a significant porosity change close to the interface. We use the planar average 
of the porosity as it is obtained in the DNS, therefore, there is no need to 
explicitly set the interface position. Since the flow is within the Stokes regime, 
the Forchheimer term in Eq. HID is neglected. 



Fig. 9. A comparison between the planar average of the stream-wise velocity obtained 
by DNS and the homogenized model, ReD — 2. 


Fig. [9] shows the results of the planar average stream-wise velocity for the 
DNS solution and the GLBM. Although the porosity, permeability, fluid prop¬ 
erties and driving forces are the same, the standard GLBM homogenized model 
over-predicts the velocity in the transition zone. The dashed line shows the ho¬ 
mogenized model that only takes the Darcy force into account. These two men¬ 
tioned homogenized models use a viscosity in the porous region which is equal 
to the free flow region. We propose to use the GLBM homogenized model but 
with a viscosity in the porous region depending on the porosity by y e g = y/e. 
As we can observe in the porous region, the latter model can perfectly predict 
the DNS result. 
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5 Conclusion 

We presented three different approaches to simulate the interaction of free flow 
with porous media flow, namely, direct pore-scale simulations, as well as ho¬ 
mogenized single-domain and two-donrains approaches. The lattice Boltzmann 
method is employed both, for obtaining the pore-scale reference solution, and 
for solving the computationally more appealing homogenized problems. 

For the two-donrain approaches, four different interface conditions for deal¬ 
ing with the physical transport through a sharp interface have been evaluated. 
Our comparison yields that the two-donrain techniques are quite sensitive to the 
interface position. To further investigate this effect, we examined two definitions 
for the interface position, i.e., the exact and the apparent position assumptions. 
However, as our results indicate, both approaches fall short with respect to ac¬ 
curacy in the vicinity of the interface if the exact interface geometry is unknown. 
As an alternative approach we consider a homogenized one-domain model that 
is based on the idea of a smooth transition zone between the free flow and porous 
media models. A simple porosity-dependent rescaling of the viscosity allows us to 
accurately reproduce the results obtained by averaging the pore-scale solution. 

In future work we aim to investigate the combination of both approaches to 
allow for the treatment of more general situations in a two-scale fashion. Since the 
discussed lattice Boltzmann schemes are suitable for REV-scale computations, 
and are also highly scalable for pore-scale simulations, they lend themselves well 
for leveraging the power of massively parallel computing architectures. 
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